R code modified from https://github.com/mjacox/Thermal_Displacement/blob/master/define_heatwaves.m
Note: “00_OISST_include_plotx.R” contains included library and gplotx function in previous Rmd files.
yrng <- seq(1982,2022)
clim_years = seq(1982,2011) #for climatology
climyrs<- clim_years[length(clim_years)] - clim_years[1] + 1 #30 yrs
monstr <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
endt <- "2022-05-22"
trackdate <- seq.Date(as.IDate(endt)-6, as.IDate(endt), by="day")
curryr <- year(as.IDate(endt))
currmo <- month(as.IDate(endt))
#### Define Marine heatwave is monthly anomaly > 90% historical anomaly
Thresh <- 0.9
Detrend <- TRUE #FALSE
if (Detrend) {
prex <- "../data_src/oisst/monthly_anom_detrend/"
outd <- "../data_src/oisst/monthly_heatwave_detrend/"
} else {
prex <- "../data_src/oisst/monthly_anom_icemask/"
outd <- "../data_src/oisst/monthly_heatwave/"
}
Recalculate_thre <- FALSE
if (Recalculate_thre) {
# It's a 3-month running window, and use quantile(probs=Thresh=0.9) to aggregate the window and get the threshold(percentile>0.9) value
sty <- future_lapply(1:12, function(j) {
if (j==1) {
winj = c(12, 1, 2)
} else if (j==12) {
winj = c(11, 12, 1)
} else {
winj = c(j-1, j , j+1)
}
monj <- fifelse(winj<10, paste0("0",winj), paste0(winj))
jstr <- monstr[j]
for (i in clim_years) {
if (winj[1] == 12 & (i-1)>=clim_years[1]) {
x0 <- read_stars(paste0(prex, i-1, monj[1], "_anom.nc"))
} else if (winj[1] == 12 & i==clim_years[1]) {
x0 <- NULL
} else {
x0 <- read_stars(paste0(prex, i, monj[1], "_anom.nc"))
}
x1 <- read_stars(paste0(prex, i, monj[2], "_anom.nc"))
if (winj[3] == 1 & (i+1)<=clim_years[length(clim_years)]) {
x2 <- read_stars(paste0(prex, i+1, monj[3], "_anom.nc"))
} else if (winj[3] == 1 & i==clim_years[length(clim_years)]) {
x2 <- NULL
} else {
x2 <- read_stars(paste0(prex, i, monj[3], "_anom.nc"))
}
if (is.null(x0)) {
x <- c(x1, x2)
datex<- c(as.Date(paste0(i, monj[2], "01"), format="%Y%m%d"),
as.Date(paste0(i, monj[3], "01"), format="%Y%m%d"))
} else if (is.null(x2)) {
x <- c(x0, x1)
datex<- c(as.Date(paste0(i, monj[1], "01"), format="%Y%m%d"),
as.Date(paste0(i, monj[2], "01"), format="%Y%m%d"))
} else {
x <- c(x0, x1, x2)
datex<- c(fifelse(winj[1]==12, as.Date(paste0(i-1, monj[1], "01"), format="%Y%m%d"), as.Date(paste0(i, monj[1], "01"), format="%Y%m%d")),
as.Date(paste0(i, monj[2], "01"), format="%Y%m%d"),
fifelse(winj[3]==1, as.Date(paste0(i+1, monj[3], "01"), format="%Y%m%d"), as.Date(paste0(i, monj[3], "01"), format="%Y%m%d")))
}
if (i == clim_years[1]) {
styx <- x
datey<- datex
} else {
styx <- c(styx, x)
datey<- c(datey, datex)
}
}
print(paste0("Before process quantile, month: ", monj, " ..."))
names(styx) <- rep(jstr, length(names(styx)))
styx <- merge(styx) %>% st_set_dimensions(3, values = as.POSIXct(datey), names = "time") %>%
aggregate(by=paste0(climyrs, " years"), FUN=quantile, na.rm=TRUE, probs=Thresh, names=FALSE)
names(styx)[1] <- jstr
styx <- styx %>% dplyr::select(jstr) %>% adrop
return (styx)
})
save(sty, file="../data_src/stats/sst_1982_2011mhw_threshold_nc.RData")
} else {
load("../data_src/stats/sst_1982_2011mhw_threshold_nc.RData")
}Recalculate_mhw <- FALSE
if (Recalculate_mhw) {
for (i in yrng) {
mmx <- fifelse(i==curryr, currmo-1L, 12L)
for (j in 1:mmx) {
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
print(paste0("Now in Year-month: ", i," - ", monj, " to get heatwave"))
z <- read_stars(paste0(prex, i, monj, "_anom.nc"))
names(z)[1] <- "heatwave"
zt <- matrix(rep(0, len=1440*720), nrow = 1440)
zt[which(z[[1]]>=sty[[j]][[1]] & z[[1]]>0)] <- 1 ## NOTE: HERE different with ref code, add criteria: ANOMALY > 0 (for MHWs)
z[[1]] <- zt
filet <- paste0(outd, i, monj, "_heatwave.nc")
write_stars(z, filet)
}
}
}ht <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/202204_heatwave.nc"))
names(ht)[1] <- "heatwave"
dht <- as.data.table(ht) %>%
.[,`:=`(x=fifelse(x>180, x-360, x))] %>% .[,.(x, y, heatwave)]
ghx <- gplotx(dht, "heatwave", returnx=TRUE, minz=-0.4, maxz=1, maxlim=180, no_coord=TRUE, legend_label="2022-04 Marine Heatwaves ", legend_pos=c(0.65, 0.09))
ghx <- ghx + geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3)
ghxzk <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "02", "_anom.nc"))
zt <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "02", "_anom.nc"))
names(zk)[1] <- "anom"
names(zt)[1] <- "anom"
zkt <- as.data.table(zk) %>% .[,`:=`(x=fifelse(x>180, x-360, x))] %>% .[,.(x, y, anom)]
ztt <- as.data.table(zt) %>% .[,`:=`(x=fifelse(x>180, x-360, x))] %>% .[,.(x, y, anom)]
#print(range(na.omit(zkt$anom))) #-9.520356 7.962329
#print(range(na.omit(ztt$anom))) #-5.414660 7.204932
gt <- gplotx(ztt, "anom", returnx=TRUE, minz=-5.5, maxz=7.5, fillopts=NA, legend_label="2020-02 Anomaly ", legend_pos=c(0.65, 0.09), maxlim=180)
gk <- gplotx(zkt, "anom", returnx=TRUE, minz=-10, maxz=7.5, fillopts=NA, legend_label="2020-02 Detrended Anomaly ", legend_pos=c(0.65, 0.09), maxlim=180)
layt <- rbind(c(1,1),c(2,2),c(3,3))
grid.arrange(gt, gk, ghx, layout_matrix=layt)if (!require("layer")) {
if (!require("remotes")) install.packages("remotes"); library(remotes)
remotes::install_github("marcosci/layer")
library(layer)
}
Reload_hw_nc <- FALSE
if (Reload_hw_nc) {
hw202002 <- as.data.table(ht) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
#.[,.(longitude, latitude, heatwave)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
ht1 <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/201002_heatwave.nc"))
names(ht1)[1] <- "heatwave"
hw201002 <- as.data.table(ht1) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
ht2 <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/200002_heatwave.nc"))
names(ht2)[1] <- "heatwave"
hw200002 <- as.data.table(ht2) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
ht3 <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/199002_heatwave.nc"))
names(ht3)[1] <- "heatwave"
hw199002 <- as.data.table(ht3) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
save(hw199002, hw200002, hw201002, hw202002, file="../data_src/stats/heatwave_1990_2020_stack.RData")
} else {
load("../data_src/stats/heatwave_1990_2020_stack.RData")
}
tl0 <- layer::tilt_map(hw199002)
tl1 <- layer::tilt_map(hw200002, y_shift=75)
tl2 <- layer::tilt_map(hw201002, y_shift=150)
tl3 <- layer::tilt_map(hw202002, y_shift=225)
map_list <- list(tl0, tl1, tl2, tl3)#palettem: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
layer::plot_tiltedmaps(map_list, palette = "cividis")